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ABSTRACT 

We study the collapse of protoclusters within a giant molecular cloud (GMC) to determine the 
conditions under which collapse is significantly disrupted. Motivated by observations of star forming 
regions which exhibit flattened cloud structures, this study considers collapsing protoclusters with disk 
geometries. The collapse of a 10 3 M Q protocluster initially a distance of 2 — 10 pc from a 10 3 — 10 6 Af Q 
point mass is numerically calculated. Simulations with zero initial relative velocity between the 
two are completed as well as simulations with relative velocities consistent with those observed in 
GMCs. The results allow us to define the conditions under which it is safe to assume protocluster 
collapse proceeds as if in isolation. For instance, wc find the collapse of a 1O 3 M0 protocluster will 
be significantly disrupted if it is within 2 — 4 pc of a 10 4 Af© point mass. Thus, the collapse of a 
10 3 M Q protocluster can be considered to proceed as if in isolation if it is more than ~ 4 pc away 
from a 10 4 M e compact object. In addition, in no portion of the sampled parameter space does the 
gravitational interaction between the protocluster disk and the massive particle significantly disperse 
the disk into the background GMC. We discuss the distribution of clusters of young stellar objects 
within the Perseus and Mon R2 star forming regions, which are consistent with the results of our 
simulations and the limitations of our results in gas dominated regions such as the Orion cloud. 
Subject headings: methods: N-body simulations, stars: formation, ISM: clouds, ISM: structure 



1. INTRODUCTION 

Over the last three decades it has become evident that 
molecular clouds serve as the site of all star formation 
within galaxies with the result that a clear understand- 
ing of molecular cloud structure and evolution is nec- 
essary to fully describe the initial conditions for star 
formation. Typical giant molecular clouds (GMCs) are 
~ 10 — 50 pc across and have masses betwe en 10 6 and 
10 7 M g (|Liszt et alJl!981l:ISolomon et al.lll987f) . They are 
nonuniform, composed of high density clumps with sizes 
of ~ 1 — 10 pc and masses of ~ 1 2 — lQ b M Pl embed- 
ded in a lower density background (|Kramer "etal]H99l 
iHever fe Terebevill998l ). These high density clumps are 
the sites of star formation within GMCs. 

Although GMCs do not globally collapse, many 
of the small-scale high density clumps undergo lo- 
cal collapse and fragmentation. The most mas- 
sive of these clumps produce clusters of young stars 
(jBallesteros-Paredes et al.l l2007t ). Throughout this pa- 
per, the term "protocluster" will refer to these mas- 
sive, high density regions within molecular clouds which 
have the potential to produce clusters of stars. The 
process by which protoclusters fragment into collapsing 
substructures and eventually protostars is complex and 
not completely understood. It may be due to a com- 
bination of Jeans instabi l ities, the decay of turbulence 
(jKlessen fc Burkertl |2000L l200l[ Kand the decoupling of 
fluid and MHD waves (jMversI Il998f) . Observations of 
star forming regions reveal that newly formed clusters of 
stars are not isolated within molecular clouds. Rather, 
star forming regions of molecular clouds tend to contain 
multiple clusters and smaller groups of young stellar ob- 



jects near each other (jAUen et al.ll2007t ). 

Since GMCs are highly nonuniform in structure and 
there is evidence of clusters forming near other clusters, 
a relevant question in cluster formation is: under what 
conditions can the collapse of a protocluster be treated 
in isolation, i.e. ignoring the influence of the surrounding 
GMC environment? 

This question is analogous to the question of how 
star and planet formation in cluster environments differs 
from formation in isolation. Studies that have consid- 
ered star and planet formation within clustered environ- 
ments have shown that environment may have modest to 
significant effects on the formation processes depending 
on the properties of cluster considered. Specifically, in 
very dense massive clusters frequent close encounters be- 
tween stars may disrupt disks and young solar systems 
and UV radiation from massive sta rs may cause pho- 
toevaporation of protostellar disk s (jBonnell fc Kroupal 
fl998t iStorzer fc Hollenbachl [1991 . However, in more 
intermediate-sized clusters, with parameters typical of 
those observed within ~ 1 kpc of the sun, interactions 
between stars are much less common, and the clus- 
ter's effects on plane t-forming disks are relatively modest 
(| Adams et al.l l2006). In this study, we consider a similar 
question, but on a larger scale and at an earlier time in 
cluster development. 

The above question can be addressed in a gross way 
with simple tidal force estimates, i.e. calculating the 
Roche limit. However the analytic Roche limit does not 
provide the details of how tidal distortion affects the evo- 
lution of the protocluster structure, how this evolution 
varies with the initial cloud structure, nor how it depends 
on the relative motion of the protocluster with respect to 
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its environment. To address these questions, one must 
complete numerical simulations of the collapsing proto- 
cluster within its GMC environment. 

To this end, we perform an ensemble of iV-body simu- 
lations to calculate the collapse of a protocluster in the 
presence of a nearby massive object in the same GMC. 
The point mass is intended to represent another proto- 
cluster or high density region that is sufficiently massive 
and compact enough to be roughly modeled as a point 
mass. These calculations treat only the gravitational in- 
teractions and ignore the effects of gas pressure, stellar 
winds, and radiation. Molecular cloud maps of star form- 
ing regions contain numerous dense gas structures with 
varying geometries, so a full numerical investigation of 
this question should consider protoclusters with spheri- 
cal, flattened, and filamentary initial geometries. In this 
study, we focus on flattened disk-shaped protoclusters. 

Our choice of protocluster geometry is partially mo- 
tivated by observations of star forming regions reveal- 
ing flattened, layered, and/or filamentary cloud struc- 
tures. For example shells of molecul ar gas surround- 
ing O B sta rs have been iden t ified by iDeharveng et al.1 
([2005ft and IChurchwell et all (|2006l ). Within many of 
these shells triggered star formation appears to be occur- 
ring (|Churchwell et al.ll2006l ; [Zavagno et al.ll2006t) . Fur- 
thermore, the distribution of young stars within embed- 
ded clusters is often aligned with the elongation or fila- 
ment structure of t he embedding gas (I Allen et alj [20071 ; 
iKumar et all 120071 : iTeixeira et all l2006h . Observations 
such as these suggest that "reduced dimensionality" may 
be a characteristic of star forming clouds. 

Flattened cloud st r uctures may be formed by OB winds 
([Weaver et all 119771: IWhitworth et all 1 1994ft or expand- 
ing H II regions (jElmegreen fc Ladalfl977ft sweeping up 
and condensing nearby cold molecular gas. Cluster form- 
ing clouds that are believed to have been compressed 
by these m echanisms include Orion B near the OBI 
association (IWilson et all 12005ft and the PR 21 ridge 
near the Cyg OB 2 association ( Schneider et alj 120061; 
IKumar et alll2007ft . lElmegreenl (|l998ft and IWhitworth! 
(|2005ft provide more detailed discussions of proposed flat- 
tening mechanisms capable of producing the structure 
observed in many sta r form ing regions. Simulations by 
iBurkert fc Hartma nn (2004) demonstrated that even in 
isolation, finite self-gravitating sheets collapse into dense 
structures with interesting geometries. 

In addition to the observational evidence for flattened 
initial states, we choose to simulate flattened protoclus- 
ters because they are harder to disrupt than spherical 
protoclusters of the same mass and radius. This is be- 
cause the disks are more centrally concentrated and thus 
more tightly bound. Therefore they are more likely to 
succeed in forming stellar clusters than their 3D coun- 
terparts, all else being equal. 

We simulate the collapse of a disk-shaped protoclus- 
ter in the presence of another dense object in the GMC 
to determine the survivability of collapsing protoclusters 
in GMCs. In § [2] the iV-body collapse calculations are 
described. The results from these simulations are sum- 
marized in § [3l and a discussion of these results in the 
context of observed star forming regions is presented in 
§ |4l The Appendix contains a discussion of the simula- 
tions of collapsing systems with analytic solutions used 
to estimate the uncertainties in A-body calculations of 



the dynamics of gaseous systems. 

2. NUMERICAL CALCULATION OF DISK COLLAPSE 

In this study, we complete an ensemble of A^-body sim- 
ulations to study the collapse of flattened protoclusters 
in the presence of a nearby dense object represented by 
a massive point particle. This dense object may be an- 
other dense protocluster, a region that is currently form- 
ing stars, a stellar cluster, or another significant den- 
sity enhancement in the GMC, provided that it is more 
massive than the protocluster and small enough that its 
gravitational influence on the protocluster can be well 
approximated as a point mass. A modified version of the 
NBO DY2 direct N-body integration code (|Aarsethll200ll . 
|2003|) is used to complete the simulations. The NBody2 
code was modified to implement a sink cell algorithm as- 
sociated with the particle initially located at the center 
of the protocluster disk. This algorithm allows particles 
that collapse to the center of the disk to accrete onto it. 
A more detailed discussion of the sink cell algorithm in 
included in the Appendix. 

A fluid dynamics code will always model a fluid system 
with more precision than will an A^-body code because 
on small scales close encounters between point masses de- 
viate from fluid behavior. For this reason A^-body tech- 
niques are not widely used to study fluid problems. Two 
significant concerns arise when modeling a fluid system 
with an A^-body code. First, A-body simulations do not 
take into account the pressure forces between point par- 
ticles. Secondly, close gravitational scattering encoun- 
ters between point particles are not physically realistic 
in fluid systems. If these two concerns are appropriately 
addressed, an A^-body code may be competitive with a 
full fluid simulation when studying the large scale gravi- 
tational evolution of a fluid. 

Specifically, if the fluid can be approximated as cold 
and pressure free, the dominant force between fluid ele- 
ments is gravity and thus the fluid elements may be mod- 
eled as individual particles moving under the influence of 
their mutual gravitational potential. Proto clusters are 
cold with typical tempera tures of < 20 K ([Liszt et al.1 
119811 : ISolomon et al.|[l987ft and can be approximated as 
pressure free fluids until the late stages of collapse. Our 
choice of radially symmetric ring structure was chosen 
specifically to minimize the number of close scattering 
interactions that particles undergo during disk evolution 
(For a more detailed discussion, see Appendix.) Thus we 
can carefully use an Af-body code to address the large 
scale question of when the collapse of a protocluster may 
be treated as if it occurs in isolation, neglecting the ef- 
fects of the surrounding GMC environment. 

In addition, for modest values of N (in this study, 
A~ < 1000) A r -body simulations are computationally 
more efficient than fluid codes. Therefore, judicious use 
of an A^-body code with a relatively modest investment 
of computational time can return a broad understanding 
of a system over a wide range of parameter space. 

To better understand the regime in which the A^-body 
simulation reliably models gaseous disk collapse behav- 
ior and thus can be used effectively in this study, the 
code must be systematically tested. Therefore, this par- 
ticular A-body code was used to complete simulations 
of collapsing systems for which analytic solutions exist. 
Specifically, the collapse of an isolated Maclaurin disk 
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and Mestel disk were simulated and the results of the 
simulations compared to the analytic collapse solutions. 
We find that the iV-body simulations predict free fall 
times for the Maclaurin and Mestel disks with errors less 
than 3.4% and 2.4% and cumulative mass profiles with 
errors less than 6.1% and 7.4%, respectively. A complete 
discussion of these test calculations and comparison to 
the pressure-free analytic collapse solutions is provided 
in the Appendix and summarized in Table [T] 

It is important to stress that this approach investi- 
gates only the gravitational interactions within the disk 
and between the disk and the point mass. Other physical 
processes in gas such as pressure and turbulence are not 
included in our simulations of protocluster collapse. The 
detailed structure of the final cluster forming cloud and 
the amount of star formation which will subsequently 
occur depends on these processes. A more in depth dis- 
cussion of how these physical processes may affect the 
results of our simulations is reserved for § 14.21 

We assume that the flattened protoclusters are formed 
when a three dimensional dense cloud is compressed into 
a two dimensional disk. This assumption is motivated 
by evidence that flattened structures in star forming re- 
gions may be created by stellar winds or shocks sweep- 
ing material into a layer. Thus the surface density of 
the protocluster disk depends on the assumed density of 
the pre-flattened three dimensional cloud from which it 
is formed. Consider a uniform density sphere which is 
compressed along one dimension. The resulting density 
distribution, a Maclaurin disk, has the surface density 
profile 
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Q<r <r d 
r > r d , 



(1) 



where r d is the disk radius (jBinnev fc Tremaineill987T ). 

The majority of our simulations considered disks that 
are initially described by Maclaurin surface density pro- 
files. These disks are weakly centrally condensed and 
consistent with the disk formation scenario described 
above. It is difficult to form a less centrally condensed 
disk distribution by flattening a realistic three dimen- 
sional clump of gas. However, flattening of a somewhat 
centrally condensed three dimensional clump of gas will 
produce an even more strongly centrally concentrated 
disk. To investigate the evolutionary differences central 
concentration may create, a limited number of strongly 
centrally concentrated disks were also completed. Specif- 
ically we considered the pre-flattened three dimensional 
dense cloud with the density distribution of a singular- 
isothermal sphere. When flattened it produces a Mestel 
disk described by the surface density distribution 



S(r) 



E(A)A 



< r < oo, 



(2) 



wh ere S(-A) is t he surface density at the fiducial radius 
A (|Mestellll963h . A comparison of the Maclaurin and 
Mestel disk evolution is included in § 14.11 

The free fall time for a particle initially at rest at radius 
r < rd in a Maclaurin disk of mass Md is 



tff(r) 




Therefore, a Maclaurin disk collapses all at once. The 
free fall time for a particle initially at rest at radius r in 
a Mestel disk is 



tff(r) 



4 y/GY,(A)A 



(4) 



(3) 



resulting in an inside-out collapse. 

Each simulated disk is composed of a central point 
mass and 25 rings of 32 equal mass particles. The mass of 
the central particle and of each ring is varied to produce 
the desired surface density distribution. Each disk has a 
total mass of M<j = 10 3 M Q and a radius of r d = 1 pc. 
(The simulated Mestel disk is truncated at rd-) For star 
formation efficiencies of ~ 0.3 these parameters are con- 
sistent with observations of young clusters with a stellar 
masses of a few hundred solar masse s and radii on the or - 
der of a parsec (|Lada fc Ladall2003t IPorras et al.ll2003h . 
For Maclaurin and Mestel disks with these parameters 
the free fall time of a particle initially at rest at rd is 
0.342 Myr and 0.525 Myr respectively. 

In this study we investigate the disruption of a proto- 
cluster disk due a massive point particle M c which rep- 
resents a neighboring high density region in the GMC a 
distance d from the center of the disk and in the plane 
defined by the disk. We choose this initial geometry be- 
cause the most efficient way to tidally disrupt an axisym- 
metric disk is by placing the massive perturber in the 
plane of the disk. A point mass coincident with the sym- 
metry axis will encourage disk collapse or, if the separa- 
tion distance is small enough, accretion onto the massive 
particle. A point mass that is neither coincident with the 
symmetry axis nor in the disk plane will tidally distort 
the disk, but not as strongly as a point mass at the same 
distance but located in the plane of the disk. 

We consider disks that have both zero and nonzero 
initial velocities with respect to the nearby point mass, 
point masses within the range M c = 1 — 10 3 Md (10 3 — 
10 6 Mq), and separation distances of d = 2 — lOr^ (2 — 
10 pc). The majority of the simulations consider disks 
with Maclaurin density profiles. A limited number of 
simulations were also completed using the more centrally 
condensed Mestel density distribution. The results for 
the Maclaurin and Mestel disks are compared in §[3] Disk 
collapse begins at time t = and continues until the disk 
is accreted onto the neighboring point mass or reaches 
a maximum mean surface density. This maximum value 
occurs at approximately the disk's free-fall time and is 
explained in more detail in § 13.31 

By modeling the nearby high density region as a point 
mass, we overestimate its effect on the collapsing disk. 
The point mass representation is only exact for density 
enhancements that have spherically symmetric mass dis- 
tributions (and a handful of special non-spherical mass 
distributions). The high density clump in the GMC will 
likely be a distributed mass which does not fall into one 
of these special cases, and will exert a force on the disk 
that is less than the force of a point particle with a mass 
equal to the mass of the clump and located at the clump's 
center. Therefore, our estimates of tidal disruption may 
not be applicable in regions containing larger amounts 
of distributed dense gas, for example the Orion (L1641) 
region (see § [4] for further discussion) . 

A large body of work exists concerning the destructive 
effects galaxy and GMC environment can have on stellar 
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clusters. iTerlevichl (|1987t) and lTheunsl (| 19911 ) completed 
iV-body simulations to investigate the effect the galactic 
tidal field and interstellar clouds have on the long-term 
evolution of initially bound stellar clusters. More re- 
cent studies have considered the disruption of a stellar 
cluster due to interactions with passing spiral arms and 
GMCs (jGieles et al.ll2006l [20071) . In addition, there have 
been many numerical studies of the disruption of stel- 
lar clusters when the original embedding gas is violently 
remove d by winds from massive stars or nearby super- 
novae (|Lada et al.lll984t iGever fc Burkertll200H lAdamsl 
120001 ). These works consider the interaction of stellar 
clusters and stellar cluster members with their environ- 
ments while our calculations consider an earlier stage in 
the life of a stellar cluster. Specifically we determine the 
effects of GMC environment on the clouds that are pos- 
sible progenitors of stellar clusters. 

3. SIMULATION RESULTS 

The initial simulations consider the behavior of a col- 
lapsing Maclaurin disk in the presence of a massive point 
particle for various mass ratios and separation distances. 
For comparison, a limited number of simulations of col- 
lapsing Mestel disks were also completed. The results of 
these calculations are summarized below. 

3.1. Interaction Outcomes 

In each simulation, the disk evolves according to one of 
three distinct behavior patterns: collapse, elongation, or 
accretion. The disk evolution depends most sensitively 
on the mass of the neighboring point mass and the initial 
distance of separation (see § 13.21 for further discussion.) 
Figure [T] provides a face on view of disks evolving via 
different scenarios. For clarity, all of the disks depicted 
in Figure [T] are initially 3.0 pc from a massive particle, 
and only the mass of the particle is varied. 

Figure [T] (a) depicts the evolution of a disk that follows 
the collapse evolution scenario. This disk collapses to a 
point-like object with little to no distortion in the disk's 
shape. A small amount of elongation along the x-axis 
does occur. This is to be expected as it is the axis which 
intersects the center of the disk and the massive point 
particle. However, the ratio of the short axis of the disk 
to the long axis of the disk ((Ay/Ax) min , see § I3.4|) re- 
mains greater than 0.5 for all disks that evolve via the 
collapse scenario. 

Elongation refers to disks that collapse while being 
tidally squeezed. These disks result in filament-like dense 
structures. The elongation scenario is further subdivided 
into two categories: weak and strong elongation. Weak 
elongation occurs when the disk collapses into a filament 
that is shorter than the initial diameter of the disk. On 
the other hand, strong elongation occurs when the disk 
is stretched into a filament that is longer than the ini- 
tial disk diameter. The latter behavior indicates that 
the force of tidal disruption is stronger than the force of 
collapse. 

Figure [T](b) displays the evolution of a Maclaurin disk 
that undergoes weak elongation. The disk collapses as 
in the collapse scenario, but at earlier times shows sig- 
nificant distortion. The differences between strong and 
weak elongation become apparent when Figure [1] (b) is 
compared to Figure [1] (c) . A disk that is strongly elon- 
gated becomes significantly tidally stretched. Tidal in- 
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Fig. 1. — Face-on view of the collapse of a Maclaurin disk 3.0 pc 
away from a high mass particle. Each point denotes the position 
of a particle within the disk. Panel a): Collapse: undistorted 
collapse of disk to a point, point mass M c = 5Mj. Panel b): Weak 
Elongation: collapse of a disk to short filament, point mass M c = 
lOMj. Panel c): Strong Elongation: stretching of disk into a long 
filament, point mass M c = Panel d): Accretion: accretion 

of disk onto a nearby massive particle, point mass M c = lOOM^. 
The location of the point mass is marked by an X where the size 
of the X indicates the mass. For each simulation, the display time 
which most clearly represents the evolution scenario was chosen. 
The times are noted in units of the analytic free fall time of the 
disk in isolation (tjj = 0.342 Myr). The size and position of the 
disk at time t = is indicated by the circle. 

teraction, rather than gravitational collapse dominates 
the disks evolution. In addition, a strongly elongated 
disk accelerates toward the system's center of mass more 
than a weakly elongated disk, which is another indica- 
tion of the strength of the force between the disk and 
the neighboring point mass. 

The evolution of a disk accreting onto the neighboring 
point mass is shown in Figure [T] (d) . The disk is sig- 
nificantly stretched out by the tidal force of the massive 
neighbor. Moreover, it travels the initial separation dis- 
tance quickly enough that the disk does not have time to 
collapse under its own self-gravity (as in the collapse or 
weak elongation scenarios) , or reach a maximum density 
as a filament (as in the strong elongation scenario). 

3.2. Dependence of Outcome on Neighbor Mass and 
Distance 

Figure [2] displays the behavior of the Maclaurin and 
Mestel disks as a function of the influencing particle 
mass, M c , and the initial separation distance, d. The 
portions of M-d space in which the systems behave dif- 
ferently are separated by lines in the figure. Systems 
that reside in the lower right portion of the plot will 
have disks that collapse without disruption or distortion 
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from the nearby clump. This is the region of M-d space 
in which protocluster disk collapse occurs as if in isola- 
tion. The band of M-d space stretching from low mass 
and low separation distance to high mass and high sepa- 
ration distance is the region in which the disk undergoes 
elongation into a filament as it collapses. Systems resid- 
ing within this regime should form filaments stretched 
out toward the massive nearby high density regions. Fi- 
nally, systems with properties in the upper left portion 
of the M-d plot are those in which the protocluster disks 
will accrete onto the nearby massive object. We should 
not expect to find individually collapsing protoclusters 
near massive dense objects with these properties. 

A constant value of initial force between the center of 
the disk and the massive particle corresponds to a slope 
of 2 in Figure [5] due to the inverse-square nature of the 
gravitational force. The dividing lines between regions 
of different system behavior have slopes slightly steeper 
than 2 due to the geometry of the interacting system. For 
smaller values of d, the disk has less distance to travel 
before it will accrete onto the nearby point mass. Fur- 
thermore, for smaller values of d the point mass will more 
effectively squeeze the disk along the y-axis and stretch 
the disk along the x-axis (where the axes are defined as in 
Figure [1]). These processes combined cause two systems 
with equal values of M c /d 2 but different initial separa- 
tions to behave differently. Specifically, the system with 
smaller g? will be more easily distorted by the nearby 
point mass, which is consistent with the results of our 
simulations. 

Calculation of the Roche limit also provides insight 
into the behavior of a collapsing disk under the influence 
of a nearby massive object. In our nomenclature, the 
Roche limit is the distance d at which an object with 
radius and mass Md will be torn apart by a point 
mass M c (|Shulll982l ). Specifically, it is the distance at 
which the tidal force from the nearby point mass is equal 
to the force that causes the disk to collapse under its 
own self-gravity. The calculation of the Roche limit de- 
pends on the density distribution of the disk and is given 
by d R = KiMd/McY^Td where k = (8/37T) 1 / 3 for the 
Maclaurin disk and k = 2 1 / 3 for the Mcstel disk. In Fig- 
ure [2] the Roche limit is denoted by a dashed line with 
slope 3 in each panel. Note that the Roche limit ap- 
proximately corresponds to the separation between the 
weak and strong elongation. This further supports the 
distinction between the regime in which the tidal disrup- 
tion dominates the system behavior and the regime in 
which self-gravitational collapse dominates the disk be- 
havior. 

The evolution of a Mestel and Maclaurin disks are 
very similar for a given disk-to-point mass ratio and ini- 
tial separation. However, the Mestel disks are slightly 
more easily disrupted. Specifically, accretion and signif- 
icant elongation occur at more modest mass ratios for a 
given initial separation, d. Comparison of the inside-out 
collapse of a Mestel disk to the uniform collapse of the 
Maclaurin disk provides insight into why this is so. The 
free fall time for the outer portion of a Mestel disk is 
about 50% longer than that of a corresponding Maclau- 
rin disk. Therefore the exterior of the Mestel disk is 
exposed to the tidal forces due to the nearby point mass 
for a longer period of time. In addition, the Mestel disk 
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Fig. 2. — Maclaurin (upper panel) and Mestel (lower panel) disk 
collapse behavior as a function of the mass of the point mass M c 
in units of disk mass and the initial separation d in units of disk 
radius. Filled circles correspond to disks that collapse to point-like 
objects with little or no distortion due to the nearby massive point 
mass, Collapse. Filled ellipses correspond to disks that collapse 
to filament-like structures under the influence of a nearby massive 
point mass, Weak Elongation. Open ellipses correspond to disks 
that are tidally stretched to long filament-like structures under 
the influence of a nearby point mass, Strong Elongation. Open 
circles correspond to disks that accrete onto the nearby point mass 
before significant collapse occurs, Accretion. Solid lines separate 
the portions of the M -d space in which systems exhibit significantly 
different collapse behaviors. The dashed line in each panel indicates 
the Roche limit. The slopes a of the solid and dashed lines are 
indicated in each panel. 

has a longer amount of time to travel toward the mas- 
sive particle before it collapses, more easily coming near 
enough for edge accretion to occur. This difference is also 
borne out by the Roche limit calculation: the power law 
coefficient k is larger for the Mestel disk which indicates 
that tidal disruption occurs at larger distances for the 
Mestel disk than for the Maclaurin disk. It is important 
to note that although the outer regions of the Mcstel disk 
are more vulnerable to distortion and accretion, the in- 
terior of a Mestel disk is extremely robust against tidal 
distortion. Even a Mestel disk that is stretched into a fil- 
ament (strong elongation) contains a massive core which 
is a remnant of the rapid collapse of the inner portion of 
the disk. 

3.3. Disk Density Enhancement 

Figure [3] presents the maximum surface density en- 
hancement factor for each simulated Maclaurin disk over 
the course of the disk's evolution or distortion. The den- 
sity enhancement factor £(t)/£o is defined as the ratio 
of the mean surface density at time t to the initial mean 
surface density. To calculate the mean surface density 
one of the outer rings in the disk is identified as the 
boundary of the disk. The average surface density is 
then the ratio of the mass contained interior to that ring 
(including the mass of the ring itself) to the area of the 
convex hull that encircles all particles interior to the cho- 
sen ring. The convex hull of a set of points is the smallest 
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Fig. 3. — Maximum density enhancement factor as a function of 
initial disk-point mass configuration in terms of M c /d? (in units 
of Md/rfy for the Maclaurin disk simulations. Plotting symbols 
correspond to behavior of disk-point mass system as in Figure [2] 

convex set that contains all of the points in the set, and 
can be informally thought of as the polygon created by 
stretchi ng a rubber band arou nd the outside of a set of 
points (jRiplev fc Rassonlll977f K For the data presented 
in Figure [3] the boundary chosen to calculate the surface 
density enhancement factor is the smallest ring contain- 
ing at least 90% of the mass (We do not choose the ring 
that encloses 100% of the mass in order to reduce "edge 
effects" discussed in the Appendix). 

As disk collapse proceeds the value of £(i)/£o in- 
creases as the particles in the disk approach the disk's 
center. An artifact of the simulation technique is that 
some disk particles pass close to the sink cell without 
entering it and continue outward on a radial path away 
from the disk's center. Thus the mean density reaches a 
maximum (at approximately the analytic free fall time 
for the disk collapse in isolation) and then begins to de- 
crease. This is not commensurate with the behavior of 
a fluid which continues to collapse until pressure forces 
slow and eventually halt the collapse. Therefore, when a 
collapsing disk reaches a maximum density we consider 
this to be the end of the collapse phase. 

While the evolution scenarios give general information 
about the fate and morphology of disrupted protoclus- 
ters, the surface density enhancement factor gives more 
detailed information about the amount of collapse that 
occurs within the protocluster. The disks that collapse 
to a point with little to no distortion attain the high- 
est values of H max /Y,o. The Mestel disks attain values of 
Yi max /T,Q which arc significantly higher than the Maclau- 
rin disks. This again is due to the fact that the Mestel 
disk collapses from the inside-out. The collapse produces 
a massive central core which is not even slightly distorted 
due to the tidal stretching/squeezing of the nearby mas- 
sive object. For larger values of M c /d 2 the maximum 
density enhancement factor drops due to the more sig- 
nificant interactions between the collapsing disk and the 
point mass. The filamentary structures created during 
these interactions have much lower surface densities than 
the point-like structures. Disks that do not collapse onto 
themselves but instead accrete directly onto the nearby 
point mass have disk surface densities that decrease as a 
function of time, and thus attain values of T, max /T,Q ~ 1 
near the beginning of the disk's evolution. 

3.4. Disk Elongation 
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Fig. 4. — Minimum axis ratio (Ay /Ax) as a function of ini- 
tial disk-point mass configuration in terms of M c /d 2 (in units of 
Md/ r d) f° r the Maclaurin disk simulations. Plotting symbols cor- 
respond to behavior of disk-point mass system as in Figure [2] 

A final way of representing the difference between the 
evolution of disks in the presence of a nearby massive par- 
ticle is to compare the amount of disk distortion which 
occurs as the disk collapses. Because the point mass 
is located on the x-axis, tidal forces cause the disks to 
collapse along the y direction faster than along the x 
direction. Therefore, as the disk evolves the ratio of 
the disks' extent along the y-axis to its extent along the 
x-axis decreases from the initial value of 1. Figure H] 
displays the minimum axis ratio (Ay/ Ax) m i n that each 
simulated Maclaurin disk acquires during collapse. As is 
evident in the plot, the larger the initial force value, the 
more effectively tidal forces squeeze and stretch the disk 
as it collapses. The scatter in this plot is due to the geo- 
metric considerations discussed in § 13. 21 which cause sys- 
tems with small separation distances to be more severely 
distorted due to tidal squeezing. Because the outer por- 
tions of a Mestel disk are bound less tightly to the center 
of mass of the disk, Mestel disks develop significantly 
smaller values of (Ay / Ai) m .j„ than do the Maclaurin 
disks with the same values of M c /d 2 . 

3.5. Disk Dispersal 

It is important to note that even a disk which is 
substantially distorted due to the interaction ultimately 
forms a dense structure. For the mass ratios and ini- 
tial separations considered in this study, the tidal forces 
are not violent enough to disperse the disk mass into 
the background GMC. Instead the tidal forces determine 
the resulting geometry of the dense structure. Even the 
most violent interactions between a disk and a point mass 
result in the disk accreting onto the massive particle. 
We assume this point mass represents a high density re- 
gion. It is likely that this region is gravitationally un- 
stable, or becomes so with the additional accreted mass, 
and therefore will collapse into a dense structure as well. 
Thus, for the considered range of initial conditions, grav- 
itational interaction alone is incapable of making dense 
structures (i.e., the disks) significantly less dense. The 
most destructive an interaction can be is to collect nearby 
gas onto another dense structure. This effectively moves 
the dense structures around but does not destroy them. 
Consequently, under the influence of gravity alone, dense 
structures are destined to remain dense structures. 

Our simulations assume that a protocluster starts in 
an initially cold state, that is there are no internal mo- 
tions. Therefore, the initial conditions force our pro- 
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toclustcrs to be gravitationally bound objects. Thus, 
the final states of the tidally distorted protoclusters are 
also gravitationally bound. In order to disperse the mass 
contained within these dense objects a mechanism other 
than gravity must be invoked, for example winds from 
young massive stars formed within or very near to the 
protocluster. 

3.6. Moving Interactions 

Observations indicate that GMCs are not quiescent 
static objects but rather they exhibit moderate amounts 
of internal motion. For GMCs with sizes of 10 — 50 pc, 
the corresponding clump to clump velocities are on the 
order of a few kilometers per second. Therefore, in ad- 
dition to the simulations described above, we complete 
simulations of the collapse of Maclaurin disks that are 
initially moving with respect to the nearby massive par- 
ticle. 

This ensemble of simulations considers systems with 
particle-to-disk mass ratios of 10 — 1000, initial separa- 
tion distances of 5 — 10 pc, and initial relative velocities 
of 1 — 10 km s _1 . The angle between the x-axis and disk's 
velocity vector varies from 15 to 45 degrees. We found 
that for these initial conditions, the disk-point mass sys- 
tems behaved in much the same way as the corresponding 
zero velocity systems. The collapse times for the disks are 
short (0.342 Myr) so that the disks do not travel a signif- 
icant distance before collapse, and therefore the non-zero 
velocities do not significantly change the behavior of the 
disk-point mass system. 

For the 60 simulations completed, over 80% of the sys- 
tems with initial non-zero relative velocities exhibited 
behavior patterns in the same category (as defined in 
§ 13. 1|) as the corresponding system with an initial rela- 
tive velocity of zero. Systems that did not exhibit the 
same behavior as their zero velocity counterparts were 
slightly more interactive. Specifically, the system with 
M c = 1000M<j, d = 7r d , and v > km s _1 accreted 
onto the nearby point mass. The disk had a positive 
initial velocity component in the direction of the point 
mass. This effectively reduced the separation distance d 
and shifted the position of the system in the behavior 
diagram (Figure [2]) to the left. 

4. DISCUSSION 

4.1. Initial Central Concentration and Tidal Disruption 

The results of our simulations imply that central con- 
centration has a modest effect on the evolution of disks 
in the presence of a massive perturber. Specifically the 
Mestcl disks, which are strongly centrally concentrated, 
have outer regions which are more easily disrupted by the 
point mass than the Maclaurin disks. For a given disk to 
point mass ratio, Mestel disks exhibit accretion onto the 
point mass and stronger elongation at larger distances 
than the Maclaurin disks, due in large part to the longer 
free fall time of the outer portions of the Mestel disk. 

While Mestel disks have exteriors that are more eas- 
ily disrupted, their interiors collapse rapidly producing 
a dense core (§ 13. 3|) . In Mestel disks that are tidally 
stretched into filaments, the result is a very dense in- 
ner core with little elongation that is surrounded by an 
elongated lower density halo of material. On the other 
hand, Maclaurin disks which are tidally stretched pro- 
duce structures that are less centrally peaked and less 



strongly elongated. Assuming stellar cluster structure 
reflects the initial protocluster structure, the ONC is a 
good example of a dense centrally concentrated cluster 
with a relat ively lower d ensity star forming region sur- 
rounding it (|O'dellll200ll ). NGC 1333 is a good example 
of a cluster which shows moderate elongation (axis ratio 
~ 2 : 1) but a roughly constant density in the central 
regions. (iGutermuth et al.ll200l . If NGC 1333 formed 
from a flattened protocluster disk, it is likely that the 
disk was less centrally condensed and more closely re- 
sembled a Maclaurin disk than a Mestel disk. 

An important point to note is that the tidal zone out- 
side of which disks can be considered to collapse as if 
in isolation is the same for both Maclaurin and Mestel 
disks. The bottom line separating the collapse and weak 
evolution scenarios in the M— d plot (Figurc[5]is the same 
for both the Maclaurin disks (top panel) and the Mestel 
disks (bottom panel). Thus, on the largest scales, disrup- 
tion of the protocluster disk is not sensitive to the specific 
surface density distribution. However, once a protoclus- 
ter is within the tidal zone, the amount of distortion the 
disk will experience does depend on the amount of cen- 
tral concentration. The more centrally condensed disks 
have exteriors that are more easily disrupted and interi- 
ors that are more robust against tidal distortion. 

4.2. Effects of Other Physical Processes 

The effects of pressure have not been included in the 
./V-body simulations and will produce some modifications 
to the results presented here. As disk collapse proceeds, 
pressure will increase and at the latest states of collapse 
may prevent the entire disk from collapsing to a concen- 
trated point. Therefore, actual protocluster disks will 
probably not have mean surface densities that are as 
high as those reached by our simulated disks. In addi- 
tion, tidal stretching of the protocluster gas will produce 
pressure enhancements along the cc-axis which may re- 
sults in less compression in the y direction, producing 
clusters with elongations less extreme than those quoted 
in §E3 

Small deviations from a perfectly smooth initial den- 
sity distribution will be magnified during collapse and 
will likely cause fragmentation and lead to star forma t ion 
(IKlessen fe Burkertl ll2000l l200l iHeitsch et al.ll200l . If 
O and B stars are formed within the protocluster, strong 
winds produced by these stars may act to disrupt the 
protocluster from the inside out, possibly halting star for- 
mation near the O and B stars and altering the shape of 
the embeddeding gas. However, these effects take place 
at a later time in cluster evolution than is considered in 
this paper. 

Our simulations considered protocluster disks com- 
posed of particles with zero initial velocity. However, 
observations of protocluster sized cloud cores reveal un- 
structured internal velocities of ~ 1 km s _1 . These ve- 
locities are less than the escape speed from the edge of 
a 10 3 M Q disk with a 1 pc radius which is ~ 3 km s _1 , 
and below the average speed required for the protoclus- 
ter disk to be in virial equilibrium ~ 2 km s _1 . Both the 

escape speed and the virial speed increase as R^ 1 as 
the disk collapses (and Rd decreases) . Collapse proceeds 
rapidly (< 0.5 Myr) once a protocluster disk becomes 
gravitationally unstable and thus it is unlikely that a 
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significant fraction of the protocluster gas parcels will 
'evaporate' due to internal motions. 

These internal motions however, will contribute to the 
internal structuring of the collapsing cloud. The angu- 
lar momentum associated with the internal motions will 
act to rearrange the density enhancements within the 
collapsing protocluster and thus perhaps give structure 
to the forming stellar cluster. This turbulence may also 
work to enhance or inhibit star formation within regions 
of the collapsing disk. However, unless the internal mo- 
tions are well ordered so that the total angular momen- 
tum of the disk is much greater than zero the angular 
momentum from the internal motions will not be capa- 
ble of supporting the protocluster against gravitational 
collapse, or changing the gross geometry of the cluster. 

4.3. Applications in Cluster Forming Regions 

With recent advances in infrared telescopes, and most 
recently, with the Spitzcr Space Telescope there has been 
a shift in our understanding of the distribution of young 
stars within molecular clouds. These higher sensitiv- 
ity observations reveal that clusters of young stars are 
not usually isolated within molecular clouds but often 
have neighboring smaller clusters which were previously 
unidentified. Examples of regions where clusters appear 
to be forming near other clusters include the Perseus 
molecula r cloud, Orion A, and the r egion surrounding 
Mon R2 (jKirk et al.ll200l lAllenlll995l : Gutermuth et al., 
in prep.) 

Our simulations allow us to estimate the regions 
around massive clusters where we should not expect 
smaller clusters to successfully form from flattened pro- 
tocluster clouds. If we consider the system from the per- 
spective of the more massive disrupting point mass, our 
calculations allow us to define a "tidal zone" around the 
dense region or massive protocluster that the point mass 
represents. Within this tidal zone smaller clusters cannot 
form without significant disruption or distortion. The 
tidal zone's size depends on the mass and radius of the 
collapsing protocluster. Specifically, the size of the tidal 
zone around a given massive object increases as the mass 
of the collapsing protocluster decreases. In addition, the 
size of the tidal zone increases as the initial radius of the 
protocluster increases. This is consistent with the Roche 

— 1/3 

calculation where a ~ M, r&. 

Recall that the high density disrupting object is treated 
as point mass in our simulations. This overestimates the 
strength and focusing of the force the disrupting object 
has on the collapsing protocluster (assuming that the 
disrupting object is not exactly spherically symmetric). 
Therefore, the estimates of the outer region of the tidal 
zones presented here are upper limits. That is, the tidal 
zones may in fact be smaller and the tidal distortion less 
apparent. 

Consider a point mass of 10 4 A/q a distance d parsccs 
away from the center of a 10 3 M© Mestel disk truncated 
at a radius of 1 pc. Figure [2] indicates that if the disk 
center lies within ~ 2 pc of the point mass, it will accrete 
onto the point mass. However, if the disk center is be- 
tween ~ 2 and ~ 4 pc from the point mass, the disk will 
collapse while being elongated into a filament-like struc- 
ture. If the disk center is more than ~ 4 pc from the 
point mass, the disk will collapse axisymmctrically, as if 



in isolation. Therefore our calculations indicate that a 
10 4 M Q cluster should not have surviving neighbor clus- 
ters of 1O 3 M within ~ 2 pc. 

The Roche limit for a 1O 3 M Mestel disk with a ra- 
dius of 1 pc in the presence of a 10 4 Mq point mass is 
d,R = 2.71pc. Therefore the Roche limit falls neatly 
inside the region where disk disruption is significant. 
However, our simulated disks exhibit some distortion at 
slightly larger distances (~ 4 pc) than the analytic calcu- 
lation. This is because the Roche calculation ignores the 
effect of the disk accelerating toward the massive particle 
as it collapses, and thus underestimates the effect of the 
gravitational force on the protocluster. 

The Maclaurin disk is even more robust against accre- 
tion onto the nearby point mass due to its shorter free 
fall time. A 10 3 Afo Maclaurin disk with radius 1 pc will 
accrete onto a nearby 10 4 M Q point mass only if the sep- 
aration distance is less than ~ 1.5 pc. If the separation 
distance is between ~ 1.5 pc and ~ 4 pc the Maclaurin 
disk will collapse into a filament-like structure, whereas 
if the separation distance is greater than ~ 4 pc, the disk 
will collapse as if in isolation. Thus, for both the Maclau- 
rin and Mestel disks, the "tidal zone" in which collapsing 
protoclusters of 10 3 M Q are compromised is ~ 4 pc from 
a 10 4 M Q high density region. 

Estimates of this kind are applicable when we con- 
sider the distribution of clusters within molecular clouds. 
For instance, these estimates are consistent with obser- 
vations of the Perseus molecular cloud which reveal the 
cluste r B5 forming in c lose proximity to the IC 348 re- 
gion. iKirk et all (|2006t ) presented an extinction map of 
the Perseus star forming region derived from the Two 
Micron All Sky Surv ey images created as a part of the 
COMPLETE survey (jRidge et al.|[2006t) . They identified 
11 "super cores" within this map, including the B5 and 
IC 348 star forming regions. These large, high density 
cores are similar to the types of protoclusters our simu- 
lations considered. 

In this region a 10 3 M Q cluster (B5) appears to be form- 
ing near the 2 x 10 3 M© IC 348 cluster. According to 
Figure [21 we should expect to find a 10 3 M Q cluster sur- 
viving near a 2 x 10 3 A#q only if it is at least 2.5 — 3 pc 
away. Assuming a distance to Perseus of 250 ±50 pc, the 
separation of IC 348 and B5 is approximately 4.5 ± 0.9 
pc away (in projection) and therefore B5 is safely outside 
of the tidal zone surrounding IC 348. 

Another region where our estimates are applicable is in 
the region containing Mon R2. Observations of this re- 
gion indicate smaller clusters surrounding the main Mon 
R2 cluster. The main Mon R2 cluster contains a total 
mass (in stars and gas) of ~ 1700Afo and has an effec- 
tive rations of 3 pc (assuming an average stellar mass 
of 0.5Mq and an extinction to mass conversion of 1 Ay 
pc~ 2 = 15M Q ). Therefore, from Figure [5] clusters with 
masses greater than 1O 3 M should not form within about 
1 — 2 pc of the main Mon R2 cluster. No clusters are 
found within 1 — 2 pc of the main Mon R2 cluster. In 
fact, the nearest region forming groups of stars is ~ 3.7 
pc (in projection) to the east of the center of Mon R2 and 
has a mass of ~ 135M Q and a radius of 1 pc (Gutermuth 
et al., in prep). Using the Roche limit calculation, wc find 
that the tidal zone in which a 135Mq cluster with radius 
may be disrupted in the presence of a 1700Mq point 
mass is 2 — 3 pc depending on the disk density distri- 
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bution chosen. Assuming the results of our simulations 
scale down in mass as the Roche calculation does, the 
estimated tidal zone is somewhat larger, ~ 4 pc. Thus 
we predict that a protocluster at a distance of 3.7 pc 
should collapse as if in isolation and not be significantly 
disrupted by the neighboring 1700M Q cluster. 

In regions containing significant amounts of distributed 
gas, such as the region in L1641 south of the Orion Neb- 
ula Cluster (ONC), collapse disruption due to nearby 
massive objects (in this case the ONC) may be less ef- 
fective than estimated by these simulations. In this re- 
gion there is a string of stellar clusters that have devel- 
oped near one a nother and near the ONC ([Allenl fl995l : 
lAllen et aLll2007t ) . Dense gas distributed around the pro- 
tocluster is likely to wash out the gravitational force felt 
by a collapsing protocluster due to a nearby massive ob- 
ject. This may allow clusters to form nearer each other 
than one would estimate solely based on the mass and 
proximity of the nearest massive neighbor. Therefore, 
care should be taken when applying the estimates pre- 
sented here to significantly gas-rich regions. 

A rough estimate of the average background density 
required to change the results of our simulations can be 
made by assuming the background gas has a uniform 
density tiq. The force on a protocluster from a 1O 4 M 
point mass a distance 4 pc from the protocluster center 
equals the force on the protocluster due to sphere of gas 
with radius 4 pc centered 4 pc from the protocluster if 
no ~ 750 cm" 3 (assuming the gas is composed of molec- 
ular hydrogen). 4 pc is the edge of the tidal zone for 
a 1O 3 M0 protocluster near a 10 4 M Q point mass. Thus 
background gas densities higher than this may begin to 
wash out the effects of nearby massive perturbers. How- 
ever, the densities required for the forces to be equal 
scale as no ~ d~ 3 . So for the background gas to have as 
much influence on the disk as a 10 Mq point mass 2 pc 
from the protocluster, a much higher average gas den- 
sity of no ~ 6000 cm -3 is required. Typical densitie s in 
GMCs range from 10 2 to 10 3 cm~ 3 (|Liszt et al.lll981l) . In 
the higher density regions, which often exhibit clustered 
star formation, d ensities may reach 5 — 15 x 10 3 cm -3 
(jKirk et al.M2006h . and so it is in these denser regions of 
GMCs that our simulations may overestimate the effects 
of nearby perturbers on protocluster evolution. 

Tidal zone estimates of the kind presented here may 
be applied to other rich star forming regions such as the 
Cyguns X region which contains many clusters as well 
as m any sites of current star formation (jSchneider et alJ 
20061 and references therein) . 

5. SUMMARY 

Our calculations of the collapse of protoclusters in 
GMCs provide the framework in which to identify regions 
where protocluster collapse can be assumed to proceed 
as if in isolation. The following conclusions can be drawn 
from the results of our simulations: 

1. Depending on the geometry and strength of the 
interaction, protocluster disks either collapse as if 
in isolation, are weakly or strongly elongated, or 
are accreted by the neighboring point mass. The 
outcomes of our simulations are presented in Figure 
[5] and are consistent with the Roche limit estimates 
of disruption. 



2. The amount of disk disruption due to a nearby mas- 
sive particle depends most sensitively on the mass 
of and distance to the particle and less sensitively 
on the distribution of material within the proto- 
cluster disk. 

3. The tidal zone defined by the separation between 
collapse and elongation evolution scenarios does 
not sensitively depend on how centrally concen- 
trated the disk is. Maclaurin (weakly centrally 
concentrated) and Mestel (strongly centrally con- 
centrated) disks have tidal zones that are roughly 
the same for a given disk to point mass ratio. 

4. The interior of disks that are initially more cen- 
trally condensed are more robust against tidal 
distortion than those that are less centrally con- 
densed. Therefore, for flattened protoclusters that 
formed from density distributions that were cen- 
trally peaked, it is likely that once the it becomes 
unstable against gravitational collapse, the inner 
portion of the disk will continue to collapse, even 
in the presence of another massive object. 

5. On account of both the short collapse times of pro- 
toclusters and the modest clump to clump veloc- 
ities observed in GMCs, the interaction outcomes 
do not depend significantly on the relative velocity 
between the protocluster and the point mass. 

6. In systems with conditions consistent with those 
observed in GMCs, gravitational interactions alone 
only produce density enhancements within the pro- 
toclusters. Thus, tidal stripping does not signifi- 
cantly disperse mass from the condensing disk. 
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6. APPENDIX 

In this Appendix we discuss the use of A-body simu- 
lations to calculate the dynamics of a collapsing gaseous 
disk. We provide the results from the simulated collapse 
of isolated Maclaurin and Mestel disks using an A-body 
code. Analytic treatments of collapsing gaseous systems 
provide generalized descriptions of collapse behavior by 
employing idealized geometries, density distributions, 
simplifying approximations, etc. However, analyses of 
this kind are limited by the very approximations that 
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make them possible. On the other hand, fluid dynam- 
ics codes are capable of simulating various initial geome- 
tries and may include more complex physical processes 
such as gas pressure, drag forces, or interactions with 
magnetic fields, etc. iV-body simulations are computa- 
tionally much more efficient than fluid dynamics codes 
and can be used with care to simul ate the behavior of 
gaseou s (or gas dominated) systems. iToomre fc Toomrel 
(|1972h is a classic example of using an TV-body simu- 
lation (specifically, a restricted three-body simulation) 
to study the dynamics of interacting galaxies. More re- 
cently, full TV-body simulations have been employed to 
study galaxy int eractions in the stellar-dynamical limit 
(|HernciuistJll993t) . 

Modeling the behavior of a gaseous disk with an TV- 
body code requires that care be taken when configuring 
the initial conditions and interpreting the results of the 
simulation. The initial configuration should minimize the 
scattering interactions that occur. This is accomplished 
by choosing an axisymmctric ring configuration where 
particles in each ring have equal mass. Specifically, the 
simulated disk contains a central point mass and 25 rings 
of 32 equal mass particles. The total mass of each ring 
is then varied to produce a a 10 3 M Q disk with radius 1 
pc a desired initial density distribution. 

The number of particles TV chosen to represent the pro- 
tocluster was determined by two competing factors: TV 
must be large enough that the particles represent the 
smoothness of the density distribution and small enough 
to minimize particle-particle scattering interactions. We 
found that a choice of TV ~ 800 produces a disk dis- 
tribution that was both smooth enough to model ac- 
curately the behavior of the isolated disks and sparse 
enough to minimize the number of scattering interac- 
tions which occurred during the collapse simulation. For 
example, a disk with TV ~ 500 particles had initial mass 
and force profiles with errors larger than the TV ~ 800 
disk. We compared simulations of disks with TV ~ 2000 
and TV ~ 800 particles. We found that the evolution of 
both disks were very similar during the first half of the 
collapse, but at later times the TV ~ 2000 disk had larger 
errors in their mass profiles and collapse times (when 
compared to the analytic solution) due to more frequent 
scattering interactions in the higher number density re- 
gions of the disk. 

Another consideration when calculating the dynam- 
ics of fluids with TV-body simulations is that the simu- 
lated behavior and conclusions drawn from that behavior 
should be consistent with the behavior of fluid systems. 
A consequence of using the TV-body simulation is that as 
the disk of particles collapses, particle-particle interac- 
tions will allow some particles to pass through the center 
of mass and continue on a radial path outward. This 
collapse and subsequent expansion is an artifact of the 
simulation method. In a fluid system, mass collapses di- 
rectly toward the center of mass until the pressure gradi- 
ent becomes too large and collapse is slowed. Therefore, 
to simulate accretion in the TV-body code, we have im- 
plemented a sink cell algorithm. 

This algorithm allows the central particle in the disk to 
inelastically "absorb" the mass of any particle that comes 
within the sphere of radius Tsink surrounding the central 
particle. Experimentation with values of this fiducial ra- 
dius resulted in the choice of r S i„k = 0.049 pc w 10, 000 



AU. We required that the sink cell radius be small enough 
that it encompasses only the central particle at t = 0, 
and large enough that it maximized the number of disk 
particles that entered the sink cell during the simulation. 
For the test calculations of isolated Maclaurin and Mes- 
tel disks, more than 90% of the disk particles collapse 
from their initial position within the disk directly into 
the sink cell. 

We model the continuous disk as a set of discrete parti- 
cles arranged in an axisymmctric ring structure. The ring 
spacing, ring mass, and central particle mass are chosen 
so that these disks have radially averaged mass and force 
profiles that are in close agreement with a Maclaurin and 
Mestel profiles. Specifically, the simulated Maclaurin and 
Mestel disks have errors of 2.4% and 3.2% in their ini- 
tial mass profiles, respectively. The initial force profiles 
have errors of 6.0% and 10.0%, respectively. These per- 
cent errors in the mass and force profiles are calculated 
as the r.m.s. error in the profile divided by M(r±/2) and 
F(r 1/2) where ri/ 2 is the initial half-mass radius of the 
Maclaurin or Mestel disk. 

One of the problems inherent in representing a fluid 
disks as set of discrete particles is that the inner and 
outermost portions of the particle disk suffer from "edge 
effects" caused by the discrete nature of the ring struc- 
ture. That is, these portions of the disks have mass and 
force profiles that differ from the edge of the analytic 
solutions. These "edge effects" are visible in the initial 
mass profiles depicted in Figure [HI and in the calculated 
free fall times in Figure \5[ We found that these effects 
do not propagate throughout the system but remain con- 
fined to their respective regions of the disk. Therefore, 
we can overcome these "edge effects" by describing the 
evolution of the disk without considering the inner and 
outermost regions in the description. 

Thus, throughout this paper the descriptions of disk 
behavior do not include the innermost or outermost ring 
of the Maclaurin disk. The region of interest for collapse 
behavior is n < r < r 2 5 for a disk with 25 rings, where 
r.i is the radius of the zth ring. The Mestel disk has the 
property that the equation of motion for a particle in the 
disk depends only o n the mass in terior to the ring con- 
taining that particle (|Mestellll963] ). This is not generally 
true for disk systems, but the Mestel disk is a special 
case where Gauss's law in disk geometry does hold. This 
allows us to embed the Mestel disk we are interested in 
simulating within a larger Mestel disk without changing 
the behavior of the smaller interior disk. Therefore, in- 
stead of ignoring the outer rings in the determination 
of disk behavior, we initially embed the Mestel disk in 
a larger disk with Ri ar ge = l-5rMesteZ, which effectively 
moves the "edge effects" outside of the region of interest 
(r < r Mestel) i and then consider the entire small Mestel 
disk except for the innermost ring (r\ < r < r Mestel)- It 
is worth noting that more than 90% of the mass within 
the Maclaurin and Mestel disks is contained within the 
regions of interest defined above and thus only a small 
fraction of the disk suffers from these "edge effects" . 

The test calculations of disk collapse considered iso- 
lated Mestel and Maclaurin disks with initial radii of 
rd = 1.0 pc and total masses of — 1OOOM . Analytic 
solutions to the pressure-free collapse of these systems 
exist and we compare our simulated results to the ana- 
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Fig. 5. — The free fall time for each particle in the simulated 
collapsing disk is plotted as a function of the initial radial posi- 
tion of the particle within the disk (diamonds). The upper (lower) 
panel displays the data from the Maclaurin (Mestel) disk simula- 
tion. Solid lines correspond to the analytic free fall time calculated 
from the collapse solution for each initial surface density profile. 

lytic solutions. For these initial conditions, the free fall 
time for a Maclaurin disk 0.342 Myr independent of ini- 
tial position within the disk. The Mestel disk has a free 
fall time for a test mass initially at r of i//(r) = 0.526r 
Myr. The resulting free fall times for the simulated disks 
are plotted in Figure [5] as a function of initial position of 
the particle within the disk. The free fall time for a par- 
ticle is defined as the time that the particle enters the 
sink cell. For particles that do not enter the sink cell, 
the free fall time is the time at which the particle passes 
closest to the disk's center. 

The "edge effects" are clearly seen in the free fall times 
of the inner and outer rings in the Maclaurin disk. The 
error in the calculated free fall times is 3.4% and 2.4% 
for the Maclaurin and Mestel disks respectively. This er- 
ror is calculated by dividing the r.m.s. error by the free 
fall time of a particle initially at the half-mass radius. 
By choosing to consider the regions of interest described 
above, the errors in the simulated behavior are reduced. 
For instance, if the entire disk was considered in the cal- 
culation of the free fall times, the errors for the Maclaurin 
and Mestel disks would be 4.6% and 2.7%, respectively. 

Figure [6] depicts the evolution of the cumulative mass 
profiles for the collapsing Maclaurin (top panel) and the 
Mestel (bottom panel) disks. The solid curves corre- 
spond to the cumulative mass profiles of the simulated 
disk at t = 0.0,0.2,0.4,0.6, and 0.8*//. The dashed 
curves correspond to the analytic solution for the cu- 
mulative mass profiles of the collapsing disk at the corre- 
sponding times. It is clear from this plot that discrepan- 
cies in the mass profile which are relics of the initial dis- 
crete particle configuration remain throughout the simu- 
lation. However, the error is contained in the outermost 
rings, and does not significantly affect the evolution of 
the bulk of the disk. At each time step, the mass profile 
agrees well with the analytic solution out to a fraction of 
the disk radius, and this fraction does not change signif- 
icantly over time. The errors in the mass profiles range 
from 2.4% -6.1% and 3.2% -7.4% for the Maclaurin and 
Mestel collapse simulations respectively for t < 0.8 *//, 
where the percent error in the mass profile is defined as 
the r.m.s. error in the profile divided by half of the total 
disk mass. 

The error in particle positions was also calculated as a 




Fig. 6. — Top panel displays cumulative mass distributions for 



the Maclaurin disk at times t = 0, 0.2, 0.4, 0.6, and 0.1 



(in or- 



der from right to left). Thicker curves correspond to later times. 
The solid curve corresponds to the simulated data while the dashed 
curve corresponds to the analytic solution for the cumulative mass 
distribution at the corresponding times. Cumulative mass distri- 
butions for the Mestel disk at the same times are displayed in the 
bottom panel. 



TABLE 1 
Simulation Errors 





Maclaurin Disk 


Mestel Disk 


Initial Force Profile, F(r)/m 
Free Fall Time, t f f 
Initial Mass Profile, M(r, t = 0) 
Mass Profiles, M(r,t < 0.8t ff ) 


6.0% 
3.4% 
2.4% 
2.4 - 6.1% 


10.% 
2.4% 
3.2% 
3.2 - 7.4% 



function of time. At a particular time and for a partic- 
ular ring of particles in the disk, we calculate the r.m.s. 
error in the positions of the particles with respect to the 
analytic collapse solution. The r.m.s. error is scaled by 
the distance that the ring has collapsed, resulting in a 
percent error in the position of a ring at a particular 
time. For the Maclaurin disk, the errors in the simulated 
ring positions range from 0.13 — 39.19% (with a mean of 
7.04%) if we consider all rings, and all times. However, 
the largest percent errors occurred in the inner-most 
rings and at times close to the collapse time for the ring. 
If instead, we consider the region of the disk between 
rings 3 and 23 and times less than 0.75*// the errors in 
ring positions range from 0.13% — 9.95% (with a mean 
of 5.01%). This calculation was also completed for the 
simulated Mestel disk. Considering all rings and times, 
the errors in the simulated particle positions ranged from 
0.14 - 16.88% (with a mean of 5.53%). If only the por- 
tion of the disk between rings 3 and 24, and times less 
than 0.75*// are considered, the errors in ring positions 
range from 0.14- 12.46% (with a mean of 5.11%). There- 
fore, the bulk of the particles in the simulated Maclaurin 
and Mestel disks have small positional errors throughout 
most of the simulation. 
These test simulations provide evidence that our N- 
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body code is indeed useful for studying systems of col- 
lapsing disks. We were able to reconstruct with our sim- 
ulations disks of discrete particles that behaved as the 
corresponding gaseous disk with the smoothed surface 
density distribution behaves. These simulations brought 
to light the "edge effects" that our discrete systems dis- 
played while at the same time revealing that those effects 
truly remained on the edge of the disks. Therefore, they 
can be consciously worked around by defining the region 
of interest as the inner portion of the disk. In addition, 



the particles in our simulation underwent very few par- 
ticle encounters during the time leading up to the disk 
collapse indicating that the choice of the axisymmetric 
ring distribution succeeded in eliminating the scattering 
effect of random particle-particle interactions. And fi- 
nally, our calculated initial force profiles, free fall times, 
mass profiles, and particle positions agreed well with an- 
alytic solutions justifying use of this simulation method 
for calculation of Maclaurin and Mestcl disk systems. 
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